

# Reminder - Load Post and Postsim functions!

##### Load Packages and Data #####
library(tidyverse)
library(haven)
library(MASS)
library(survey)
VSG <- read_stata("VSG_Coded.dta")

ANESGSS_Design <- svydesign(id=~1, 
                         weights=~weight, 
                         nest=TRUE, 
                         data=ANESGSS2020)
VSG_Design <- svydesign(id=~1, 
                            weights=~weight_genpop_2020Sep, 
                            nest=TRUE, 
                            data=VSG)

##### COVID Concern #####

mod_concern_VSG <- svyglm(covidconcern ~ auth_2016 + age + male + black + hispanic + education + income + unemployed + south + engagement + auth_2016*engagement, design = VSG_Design)
concern_VSG <- post(mod_concern_VSG, x1name =  "auth_2016", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)

plot(x = c(0,.25,.5,.75,1), 
     y=c(concern_VSG@est[3,1,1], concern_VSG@est[3,1,2], concern_VSG@est[3,1,3], concern_VSG@est[3,1,4], concern_VSG@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Political Engagement", 
     axes=FALSE, 
     ylim=c(-.49,.45), 
     xlim=c(0,1))
axis(2, at=c(-.6,-.4,-.2,0,.2,.4, .6), labels=c("","-.40","-.20",".00",".20",".40", ""), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(concern_VSG@est[3,2,1], concern_VSG@est[3,2,2], concern_VSG@est[3,2,3], concern_VSG@est[3,2,4], concern_VSG@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(concern_VSG@est[3,3,1], concern_VSG@est[3,3,2], concern_VSG@est[3,3,3], concern_VSG@est[3,3,4], concern_VSG@est[3,3,5]))
abline(h=0, lty=3)
title("VSG 2020", line = -.1)

##### Health Restrictions #####

mod_covidrestrictions_VSG <- svyglm(covidrestrictions ~ auth_2016 + age + male + black + hispanic + education + income + unemployed + south + engagement + auth_2016*engagement, design = VSG_Design)
covidrestrictions_VSG <- post(mod_covidrestrictions_VSG, x1name =  "auth_2016", x1vals = c(0,1), x2name = "engagement", x2vals = c(0, .25, .5, .75, 1), n.sims = 1000, seed = 1996, digits = 3)

plot(x = c(0,.25,.5,.75,1), 
     y=c(covidrestrictions_VSG@est[3,1,1], covidrestrictions_VSG@est[3,1,2], covidrestrictions_VSG@est[3,1,3], covidrestrictions_VSG@est[3,1,4], covidrestrictions_VSG@est[3,1,5]), 
     type="p", 
     ylab="Effect", 
     xlab="Political Engagement", 
     axes=FALSE, 
     ylim=c(-.49,.45), 
     xlim=c(0,1))
axis(2, at=c(-.6,-.4,-.2,0,.2,.4, .6), labels=c("","-.40","-.20",".00",".20",".40", ""), las = 2)
axis(1, at=c(0,.5,1), labels=c("Low","Medium","High"))
segments(x0 = c(0, .25, .5, .75, 1), y0 = c(covidrestrictions_VSG@est[3,2,1], covidrestrictions_VSG@est[3,2,2], covidrestrictions_VSG@est[3,2,3], covidrestrictions_VSG@est[3,2,4], covidrestrictions_VSG@est[3,2,5]), x1 = c(0, .25, .5, .75, 1), y1 = c(covidrestrictions_VSG@est[3,3,1], covidrestrictions_VSG@est[3,3,2], covidrestrictions_VSG@est[3,3,3], covidrestrictions_VSG@est[3,3,4], covidrestrictions_VSG@est[3,3,5]))
abline(h=0, lty=3)
title("VSG 2020", line = -.1)



